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A geometric method, namely the Ricci curvature criterion, is used 
to study the relative instability of evolving A'"-body gravitating sys- 
tems with massive centres. The performed numerical experiments 
show that the Ricci curvature decreases with the increase in central 
concentration, thus indicating an increasing instability for systems 
with massive centre — either in the form of a point mass or a core 
formed during the evolution of single and multimass systems. Sim- 
ulations of such systems do indeed suggest that they have faster 
gravothermal evolution than less concentrated ones. 
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1 Introduction 

It has been known for sometime [1,2] that gravitational systems display dy- 
namical behaviour that is closer to that of hyperbolic systems than to near 
integrable ones. In particular, it was shown that spherical systems display neg- 
ative configuration space curvature for the majority of initial conditions. As is 
well known [3] , the behaviour of hyperbolic dynamical systems is very different 
from that of near integrable ones. It follows therefore that the study of the 
stability properties of iV-body gravitating systems may throw light on basic 
dynamical properties of star clusters and galaxies which may not be recovered 
by standard techniques of galactic dynamics, which assume that large spheri- 
cal A^-body gravitational systems can be treated as integrable (perhaps with 
the inclusion of additive noise to model discreteness effects — even this how- 
ever is usually neglected for systems with large standard two body relaxation 
times) [4]. 
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A first step in the study of cliaotic befiaviour in numerical models is obviously 
to find proper descriptors of such behaviour and to determine what they tell 
us about a given system. The difficulties of the interpretation of computer 
results concerning A'^-body models are known since the pioneering numerical 
work of Richard Miller [5] . The rigorous study of statistical properties of many 
dimensional non-linear systems, such as the A^-body gravitating systems, is 
also associated with difficulties concerning the application of certain criteria 
and methods. For example, principal difficulties are associated with the use of 
Lyapunov exponents, given the fact that they are discontinuous functions of 
bifurcation parameters [6], and due to their discontinuity with respect to the 
computer-created images of A^-body systems [7]. Lyapunov exponents are also 
asymptotic quantities which are not particularly suitable for the study for the 
local (in time) chaotic behaviour, of A^-body gravitational systems. 

In the present study of the stability properties of the phase space trajectories 
of certain A^-body configurations we use the geometric approach, well known in 
the theory of dynamical systems [8]. In physical problems these methods have 
been introduced originally by Krylov [9]. In stellar dynamics they have been 
first used by Gurzadyan & Savvidy [1,2] (for further reviews see [10]). In the 
aforementioned work the sign indefiniteness of configuration space curvature of 
A^-body gravitating systems was also shown — meaning that, in general, these 
systems are not, strictly speaking, Anosov systems. The latter fact dictates the 
necessity of searching for more weak but checkable criteria in characterising 
the statistical properties of A^-body trajectories. One possible way is to find 
criteria of relative chaos, i.e. conditions which should describe the instability 
properties of a given system with respect to another compatible system. Such 
a criterion appears to be one involving the so-called Ricci curvature, which we 
will use in the present study. 

The question we are interested in here is the effect of the regular field produced 
by a central mass concentration on the stability properties of an evolving spher- 
ical N-body system. The massive centre in real star clusters or galaxies can 
appear as a massive black hole or central dense stellar core (resulting from 
core collapse or mass segregation resulting from energy equipartition in mul- 
timass systems). However, various physical processes directly associated with 
both these phenomena (loss cone effects, for example) and which are crucial 
in the vicinity of the centre, are out of the scope of our analysis. 

Motivation for our study stems in part from a series of results where a central 
mass concentration has been found to affect the dynamics of gravitational sys- 
tems in an important way. The destabilizing role of a regular central field has 
been observed by Gurzadyan & Kocharyan [13] using the Ricci curvature as 
a chaos indicator; in that study only small static systems were investigated, 
while in the present paper we will investigate systems with far larger number of 
particles which start near detailed dynamical equilibria and which may evolve 
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in time. Even before it has been noticed in A^-body simulations that asymmet- 
ric systems quickly evolve towards more isotropic shapes when a significant 
central mass is present ([11]; see also [12]) an efi:ect due to the destabiliza- 
tion of the main periodic orbits. The increase in the rate of relaxation for 
systems with massive centre (black hole) has been studied also by Ranch & 
Tremaine [14] where the effect of resonant relaxation was investigated. It was 
also shown that the central mass makes spherical system more unstable with 
respect to the gravothermal catastrophe [15]; that result has been obtained 
by means of methods of catastrophe theory. The stability of systems with a 
central point mass has been studied also by [16]; it was shown that not only 
the neighbouring particles have role in the stability, but also the mean field. 

The content of this paper is as follows. In Section 2 we describe the geometric 
criterion for intstability used in this paper, in Section 3 we discuss the numer- 
ical method. Section 4 presents the results which include simulations of single 
mass systems without initial central concentrations (Section 4.1) and systems 
with initial central concentrations or with different particle masses (Section 
4.2). The results are discussed in Section 5. 



2 Ricci curvature criterion 



For an introduction to the basic ideas behind the use of geometric methods 
in the theory of dynamical systems the reader is referred to reference [8]. 
Here we briefiy present the basic concepts needed in the interpretation of 
the numerical results. The Ricci curvature as a criterion for relative chaos 
has been introduced by Gurzadyan and Kocharyan in [17], where a kind of 
classification of various stellar configurations by their degree of chaos had been 
also achieved. The criterion is based on the computation of the Ricci curvature 
in the direction of the geodesic velocity u of the phase space of the system, 
and is defined as follows [19,18]: 

II " II M=i 

where Ric is the Ricci tensor, s is the geodesic time, k^,^ is the two-dimensional 
curvature, and the deviation vectors are chosen as n^-Lu and n^-Lrii, for all 

This form of the Ricci curvature arises while averaging the Jacobi equation — 
which links the dynamics of the flow with the geometry of the phase space. 
One then obtains 

^ = ~r„ + {||V„nn. 
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Then, the criterion of relative instability reads: among two systems with r\ 
and r2 within an interval of affine parameter s - [0, s*], the system with smaller 
negative r should typically be more unstable, i.e. unstable with higher prob- 
ability, where 

' = ^oMsM'^^^ ri<0;ri<r2; 

Certain loss of information on configurations of the system is obvious, however 
the Ricci curvature does contain information on typical systems, as distinct 
to the scalar curvature which does not contain it. 

The explicit form of the Ricci tensor readily follows from the Riemann tensor 
is given by [1] 

Ric^, = -{1/2W) [AWgx, + (3iV - 2)Wxp] 
+{3/AW^) [{3N - 2)WxWp] - 

where 

TV 

W^E- V{q) ^E-Y, Gmam,/rab, 

b<a 

and rrii are the masses of the particles of the A^-body system. 



3 Numerical method and choice of system 

The numerical A^-body routine used in this study is the NB0DY2 code de- 
veloped by Sverre Aarseth [20] and kindly made available by him. This is 
a direct summation code which uses individual time-steps for each particle 
in the simulation and speeds up the force calculation by using a neighbour 
scheme which separates the force calculations for neighbouring particles from 
those further off. These improvements take into account the very different 
natural times (~ l/y^, p being the local density) in a gravitational system 
and the fact that, at a given point, the irregular force due to nearby neigh- 
bours varies much faster than the regular force due to particles further off. For 
highly inhomogeneous centrally concentrated systems these refinements can 
lead to a considerable increase in efficiency while essentially maintaining the 
high accuracy of direct summation method. The errors in the calculations (as 
measured by energy conservation) are controlled by an accuracy parameter 
T] which determines the size of the integration time-steps. These errors are 
constant for values of rj below 0.01 and increase as rj'^. We have found that 
a value of rjirr = 0.02 (controlling the irregular time-step) and rjreg = 0.04 
(for the regular time-step) gave reasonable results while maintaining efficient 
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running of the code. For < 1000 the maximum number of particles allowed 
in the neighbour sphere was 10 + s/N while for larger N it was taken to be 

VV '^irr 4 ) 

The systems are taken to start from Plummer sphere configurations with the 
initial conditions determined by the method described in [21]. Using the units 
of Heggie & Mathieu [22] (by setting the total mass and the gravitational 
constant to unity while keeping the total energy of the system dX E = —0.25) 
the Plummer model is defined completely: the virial radius is fixed to unity 
and the mean virial equilibrium crossing time is equal to time units. 

The NB0DY2 code is not designed to work with zero softening [20], therefore 
a small softening, again of Plummer form, has been added to avoid numerical 
overflow due to close encounters and large angle scatterings. The softening 
parameter was taken as e — 2/N which scales as the ratio of the minimum to 
maximum impact parameters of standard relaxation theory [23]. 

Since the Ricci curvature determines the average or typical instability response 
to random perturbation, it requires not too large dispersions. It is therefore 
necessary to remove the larger contributions due to close encounters. In [24] 
the Ricci curvature time series corresponding to a given simulation have been 
averaged over small time intervals while filtering out terms that were much 
larger than the typical value. Here instead we simply not include contributions 
to the Ricci curvature due to interactions of particles within a (softened) radius 
of 0.05 units from each other (5% of the virial radius). The typical fraction of 
operations excluded by this procedure is about a few in ten thousand. 



4 Results 

4-1 Variation of the Ricci curvature as a result of core contraction 

The diagrams on the left hand side of Fig. 1 show the time evolution of the 
Ricci curvature, as a function of the crossing time, for systems starting from 
different A'^-body realizations of the same Plummer model. The uppermost 
plot corresponds to the evolution of 200 particle system, while the follow- 
ing plots correspond to the evolution of systems consisting of 600, 1000 and 
1400 particles respectively. Clearly the Ricci curvature is always negative and 
becomes more so as the systems evolve. This can be understood readily. 

Since the pioneering work of Antonov [25] and Lynden-Bell and Wood [26], it 
is well known that gravitating systems in virial equilibrium have no entropy 
maxima. Instead they can continually evolve towards higher entropy states 
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characterised by the appearance of progressively tighter cores surrounded by 
a diffuse halo. 

It is therefore natural to take the evolution of the core radius as an indicator of 
the evolutionary state of a spherical system. Plummer spheres do indeed tend 
towards more concentrated configurations as they evolve, we must therefore 
assume that these higher entropy states arc characterised by smaller values of 
the Ricci curvature, i.e. with higher dynamical entropy states. Note that here 
a correlation does exist between the dynamic entropy given by Ricci curvature 
and the thermodynamic one, which is not the case for any physical states of 
gravitating configurations because of their non-compact phase space. 

NB0DY2 contains a routine for finding the core radius of a spherical system 
using the prescription of Casertano and Hut [27] . Particles inside a given radius 
from the origin (usually the half mass radius) are selected. For each of these 
particles a list is made of the six nearest particles. The "mass density" for 
particle i is then defined as pi = 3M5/(47rr|) where rg is the distance to the 
particle farthest away from particle i in the list and M5 is the sum of the 
masses of the five other particles. The density centre of the system 



is then calculated. With the help of this quantity the core radius is then defined 
to be 



Rr 



E 



(2) 



The plots on the right hand side of Fig. 1 show the evolution of the core radius 
as calculated from Eq. (2) for the runs which correspond to the Ricci curvature 
plots adjacent to them. These diagrams clearly show that the evolution of the 
Ricci curvature properly correlates with the evolution of the core radius, in 
the sense that they both decrease together, thus confirming the conjecture 
outlined above. Only the scaling in the Ricci curvature undergoes changes 
with the number of particles, but not its behaviour. 



4.2 Massive centre and multimass systems 



As we have seen above, as the Plummer sphere models become more concen- 
trated their dynamical entropy, as measured by the negative Ricci curvature, 
increases. Hence if there is any correlation between this quantity and the rate 
of evolution of a system towards a tighter core, then we expect centrally con- 
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Fig. 1. Evolution of the core radius and Ricci curvature for different realizations of 
the same Plummer model 



centrated systems to evolve faster, thus outlining the self amplifying nature 
of gravothermal instability. 



To see if the above argument is valid we have conducted simulations in which 
one of the particles, which is initially placed at the centre of the system, has 
a mass that is much larger than the other remaining 999 particles. The mass 
of the central particle was taken to be either ten, fifty or a hundred times the 
mass of the other particles in the simulation. Other system parameters like 
the total mass, energy, size etc were kept constant. The Plummer model is 
rescaled so as to start from virial equilibrium including the heavier particle. 
The systems stay near virial equilibrium throughout the simulations — in 
the case when the central particle is ten times more massive than the other 
particles, the virial ratio T/W does not vary by more than 2%. This suggests 
that our models start from near a detailed dynamical equilibrium and that 
initial departure from such an equilibrium would not play a drastic role in 
determining the subsequent evolution. Because of its significantly larger mass 
the particle remains essentially within the inner 0.3% of the (tidal) radius of 
the system even for the case when it is only ten times the mass of the other 
particles. 
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Fig. 2. Evolution of core radius and Ricci curvature for a = 1000 system in which 
of one of the particles, initially located at the centre, is ten times more massive than 
each of the remaining 999 particles 
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Table 1 

Time-scales (in crossing times) derived from the values of scalar curvature and 
averaged over ten outputs during the first crossing time for N = 1000 systems in 
which one of the particles, originally located at the centre, has a mass Mp that 
cither 10, 50 or a 100 times more massive than the mass of each of the remaining 
999 particles of mass M. The third column estimates the RMS error in the ten 
outputs 

We show the corresponding plots for the Ricci curvature and core radius time 
evolutions in Fig. 2, which should be compared with the N = 1000 plots of 
Fig. 1. It is evident that the evolution of the core radius is much faster for the 
system with a central massive particle confirming that such system is more 
unstable. Thus, again the evolution of the Ricci curvature correlates with the 
evolution towards a denser core. This result is perhaps not that surprising since 
core contraction is related to the negative specific heat of open gravitational 
systems — which means that a subsystem losing energy increases its kinetic 
energy. This process is amplified when the velocity gradient in the mean field 
is steep, which is precisely what happens when there is a central density cusp 
or concentration. 
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Similar behaviour was found for the system with higher masses for the central 
particles. However, for N = 1000, the data for these systems were found to 
be more noisy. Both the Ricci curvature and the core radius calculations were 
found to be dominated by the position of the central massive body and by 
close encounters of this particle with the surrounding ones. Also, due to the 
continuing tightening of the core and the formation of hard binaries, numeri- 
cal difficulties were encoimtered since the NB0DY2 contains no provisions for 
handling such effects. The simulations had to be stopped because they had 
slowed down dramatically after about twenty two crossing times. It is inter- 
esting to note that, in the case when the central particle has a mass equal 
to a hundred times that of the other particles, the core radius had already 
contracted to 0.04 units after a few (~ 3) crossing times. Table 1 shows the 
characteristic time-scales (in units of the crossing time) averaged over the first 
ten outputs (corresponding to an interval of about a crossing time) for the dif- 
ferent models. These are derived from the values of the scalar curvature (which 
are less noisy). The use of scalar curvature is motivated by the analytical proof 
in references [1,2] that the relaxation time-scale of spherical A^-body systems 
is determined by the scalar curvature. 

Fig. 3 shows the behaviour of the Ricci curvature and core radii for systems 
consisting of two thousand particles when the central particle is ten, fifty and 
one hundred times more massive than the remaining particles. The estimate 
of the core radius was again found to be very noisy if the central particle is 
included. Moreover this estimate then reflected the large mass concentrated in 
the central particle rather than the overall contraction of the core due to the 
dynamical evolution. For these reasons, the contribution of this particle was 
removed when calculating the aforementioned quantity. This type of filtering 
was not applied in the calculation of the Ricci curvature (which is a charac- 
teristic of the self consistent dynamics). As a result the time evolution of this 
quantity contains larger fluctuations. Again, the larger the central mass the 
faster the evolution. While the presence of a central particle only ten times 
more massive than the remaining ones (a perturbation of 0.5%) was not found 
to have a marked effect on the evolution, systems with central masses of 1/40 
and 1/20 the total mass of the system evolved signiflcantly more rapidly. These 
results were also conflrmed for systems with A^ = 3500. The results for these 
systems are shown in Fig. 4 Once more we see that when the mass of the 
central particle is sufficiently small, the evolution is relatively slow. It appears 
that (at least for the A^ = 1000, 2000 and 3500) a central mass of ~ 1% is 
sufficient to significantly accellerate the evolution of spherical gravitational 
systems. 

In a system consisting of particles with different masses, interactions between 
particles leading to equipartition of energy means that more massive particles 
tend to lose kinetic energy and reside near the centre. This leads to a growth 
in the central mass concentration and the processes described above will then 
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Fig. 3. Evolution of core radius and Ricci curvature for N = 2000 systems in which 
of one of the particles, initially located at the centre, is (from top to bottom) is 10, 
50 and 100 times more massive than each of the remaining 1999 particles 
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Fig. 4. Evolution of core radius and Ricci curvature for N = 3500 systems when 
one of the particles, initially located at the centre, is 10 (top) or 150 hundred times 
more massive than each of the remaining 3499 particles 

take place. As can be seen from Fig. 5 this is clearly the case. Actually the 
evolution in the multimass system is even much more rapid. This is due to the 
fact that here there are two effects — the exchange of energy leading to high 
central mass and the effect of the central mass itself — which reinforce each 
other and drive the evolution. After about thirty crossing times many particles 
have effectively escaped the system and many more are ejected from the central 
areas and, except for the tightening central core, the system becomes more 
diffuse. This is clear from the evolution of the tidal radius shown in Fig. 6. 
The core radius becomes very small and only a few particles are enclosed by 
it (rendering the calculation of that quantity to be itself inaccurate). Hard 
binaries start forming and become tighter as the system evolves. These effects 
cause the integration to practically stop at about forty five crossing times. 
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Fig. 5. Evolution of core radius and Ricci curvature for N = 1000 multimass system 
with Salpeter mass function and in which the highest mass particle is ten times more 
massive than the lowest mass one 




o 1 o 20 30 4-0 



Fig. 6. Evolution of tidal radius of the system considered in Fig. 5 

Before that, the region surrounding the core reexpands as the negative specific 
heat is reversed by the energy stored in the binaries. This state of the system 
is actually more regular than the intermediate core collapse since many of 
the escaping particles are moving in regular orbits and the outer areas are 
very diffuse. This is one of the reasons for the increase of the Ricci curvature 
near the end of the simulation — the others being that at these later stages 
a large contribution to the Ricci curvature is due to the tight central binaries 
performing regular two body motion. At this stage we have reached the limits 
of the applicability of our geometric chaos indicator — an effect due to the 
non-compactness of the phase space. 

Multimass systems with N — 2000 and N — 3500 were also found to evolve 
faster than their single mass counterparts with similar global properties. This 
can be seen from Fig. 7 which shows the the evolution of Ricci curvatures 
and core radii for systems with Salpeter mass function and where the most 
massive particle has a mass ten (for the case when N = 2000) and twenty (for 
N — 3500) times larger than the least massive particle (respectively). Again 
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Fig. 7. Evolution of core radius and Ricci curvature for N = 2000 (top) and 
N = 3500 multimass systems with Salpeter mass functions. The most massive 
particles are 10 and 20 times more massive than the least massive ones in the 
N = 2000 and the N = 3500 respectively 

there is a clear general correlation between the evolution of the two quantities. 
However, as can be seen, there is some difference in the detailed behaviour. 
This might be expected since the Ricci curvature and the core radius are 
defined in rather different ways. Therefore fiuctuations due to the finiteness 
of the system are expected to be different. In particular, the core radius is 
not a characteristic quantity of the dynamics of the system but is a (rather 
arbitrarily defined) indicator of its central concentration and depends only on 
the positions of the particles. 
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5 Discussion and conclusions 



In the present study the statistical properties of certain A''-body systems have 
been analysed, particularly the role of the field of a massive centre is revealed. 
The massive centre in stellar systems can be imagined either as a central point 
mass (massive black hole) or dense stellar core formed during the gravothermal 
evolution of the system. Therefore we have considered both cases, even though 
one can expect a priori that their resulting dynamical effect on the system 
should be similar (and indeed this turns out to be the case). 

The main result of the present study is that a massive centre leads to an in- 
crease in the instability (i.e. mixing) properties of evolving N-body gravitating 
system. This increase in the mixing properties leads to faster gravothermal 
evolution for centrally concentrated systems underlining the self amplifying 
nature of the gravothermal catastrophe, in which central concentration leads 
to an increase in the temperature gradient which in turn leads to an accel- 
eration of the process. Indeed, it was found that systems with central mass 
concentration, either in the form of a point mass initially placed at the centre 
or multimass systems which evolved a central concentration because of the 
accumulation of massive particles near the centre, evolved faster. 

One prediction of this study, therefore, is that stellar clusters with initial 
central mass concentrations (such as primordial black holes for example) will 
evolve faster towards core collapse. The faster mixing rate also means that 
a perturbed system returns to its unperturbed state at a faster rate. This 
indicates that, for example, a star cluster with massive centre being shocked 
by Galactic tidal field, will recover its spherical shape faster than a cluster 
without a central mass. This problem is associated with that of the ellipticity 
of Galactic globular clusters, attracting now much attention. The relative rate 
of evolution should be proportional to the values of the Ricci curvature, which 
by definition describe not the absolute but relative degree of instability of 
different systems. 

Since centrally concentrated states of gravitating systems have higher thermo- 
dynamic entropy, it follows that the dynamical entropy (which increases with 
decreasing curvature) is correlated with the thermodynamic entropy in grav- 
itational systems (at least for the cases studied). In previous studies [24,28] 
it was shown that the existence of large scale ordered motion also decreased 
the dynamical entropy — thus also suggesting the above conclusion. In par- 
ticular, El-Zant [24] has found that the dynamical entropy increased during 
the evolution of initially ordered systems with macroscopic (e.g., plasma type 
collective) instabilities. It is important to note however that this correllation 
between thermodynamical and dynamamical entropy does not hold in the 
later stages of evolution of A^-body systems, when effects due to the non- 
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compactedness of the phase space are prominent. This is because, striclty 
speaking, the methods we have employed here are apphcable for systems with 
a compact phase space. The dynamics in these later stages are actually more 
regular as the system disintegrates. The peculiar thermodynamic properties 
of gravitational interaction allows such states to have higher entropy. 

The correlation of the exponential instability with the various physical pa- 
rameters of A^-body systems suggests that the mixing properties predicted 
by the existence of this instability may play an important role in the evo- 
lution of gravitational systems. Of particular importance is its effect on the 
time-scales of relaxation phenomena as originally pointed out by Gurzadyan 
& Savvidy [1,2]. 

One has, however, to note the following. The characteristic time-scales, which 
are found to be a fraction of a crossing time (Table 1), need not be simply 
equal to the actual time scales of relaxation and physical evolution of a grav- 
itational system. The latter properties are determined by the rate of diffusion 
in the action variables which may be much slower than the g time [29]. How- 
ever, since no universal relaxation time scale does exist for all N-body systems 
and different physical quantities may have different relaxation time scales, this 
does not imply that the exponential divergence has no physical observable ef- 
fect. Neither does it mean that the diffusion rate is determined for all action 
variables by classical relaxation theory since this assumes an integrable system 
which remains so under perturbations due to discreteness — and this contra- 
dicts the observed exponential divergence. Moreover, though the divergence 
is local, it takes place at almost all initial conditions so that the structure of 
the phase space is such that, for a given trajectory most nearby trajectories 
diverge almost everywhere. In this sense the divergence, although locally char- 
acterised, is a global phenomenon and must have an effect on the statistical 
properties of the dynamical system (see e.g. reference [29]). The criticism by 
Goodman, Heggie and Hut [30] of the linear theory of stability (Jacobi equa- 
tion) and its relation to macroscopic relaxation, is therefore removed due to 
these results from the theory of dynamical systems. Cerruti-Sola and Pettini 
[31] suggest to add the first and second order time derivatives of the potential 
energy in the Jacobi equation, which perhaps can be of some interest for some 
systems, but not for stationary stellar systems where the potential energy does 
not undergo significant time fluctuations. Also, they state that due to Schur's 
theorem, if the curvatures are constant, one can replace the two-dimensional 
curvatures with the averaged (Ricci or scalar) curvatures without loss of infor- 
mation. This, however, does not follow from that theorem. Finally, their use 
of the Eisenhart metric is also not justified, since that metric is not positively 
defined. 

These examples once more show the difficulty of the problem of characterising 
and interpreting the instability of nonlinear many dimensional system like 
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A^-body gravitating systems. Among recent results in this area one can also 
mention the indications for the possible existence of a universal scaling in 
iV-body dynamics, either of fractal nature [32] or due to Feigenbaum period 
doubling bifurcations [33]. 
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